Excitation spectra of disordered dimer magnets near quantum criticality 
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For coupled-dimer magnets with quendied disorder, we introduce a generalization of the bond- 
operator method, appropriate to describe both singlet and magnetically ordered phases. This allows 
for a numerical calculation of the magnetic excitations at all energies across the phase diagram, 
including the strongly inhomogeneous Griffiths regime near quantum criticality. We apply the 
method to the bilayer Heisenberg model with bond randomness and characterize both the broadening 
of excitations and the transfer of spectral weight induced by disorder. Inside the antiferromagnetic 
phase this model features the remarkable combination of sharp magnetic Bragg peaks and broad 
magnons, the latter arising from the tendency to localization of low-energy excitations. 



Magnetic quantum phase transitions (QPT) have at- 
tracted enormous interest over the past two decades, with 
intriguing aspects such as Bose-Einstein condensation of 
magnons, exotic criticahty, non-Fermi liquid behavior, 
and secondary instabilities, e.g., towards superconduc- 
tivity [lH3. The influence of quenched disorder, being 
inevitable in condensed-matter systems, on QPT is a less 
explored field, although a number of important theoret- 
ical results are available 0, Q : Disorder can modify the 
critical behavior or even destroy the QPT; in addition it 
can produce singular response in off-critical systems via 
quantum Griffiths physics dominated by rare regions . 

While the thermodynamics of disordered model sys- 
tems near quantum criticality has been studied using 
a variety of theoretical tools [11, rather little is known 
about the dynamics of excitations in this fascinating 
regime, mainly because numerical methods are either 
restricted to the ground state or work in imaginary 
time where real-frequency spectra are difficult to extract, 
while analytical methods are restricted to low energy and 
to the limits of either weak or strong disorder. For quan- 
tum magnets, this issue is pressing, as high-resolution in- 
elastic neutron scattering (INS) experiments are getting 
access to magnetic excitations near QPTs and 



suitable materials with disorder created via intentional 
doping are available 12- 15^. 

The aim of this paper is to close this gap: For coupled- 
dimer magnets, being paradigmatic model systems for 
magnetic QPT [iHsf , we propose a suitable generalization 
of the bond-operator method [l6[ to cases with quenched 
disorder. This enables the numerical calculation of mag- 
netic excitation spectra in both paramagnetic and mag- 
netically ordered phases, including the vicinity of the 
QPT, with disorder being treated exactly in finite-size 
systems. We apply the method to the square-lattice bi- 
layer Heis enberg model, studied in detail in the disorder- 
free case 
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-20[, with different types of exchange ran- 
domness. Near quantum criticality, we generically find a 
strong broadening in both energy and momentum of the 
low-energy magnetic response, while that at high energy 
is less affected by disorder. An interesting dichotomy 
follows upon slightly moving into the antiferromagnetic 



(AF) phase: the strongly inhomogeneous AF features 
low-energy magnons that are distinctly broadened due to 
disorder-driven mode localization, but at the same time 
displays sharp Bragg peaks in its elastic response. We 
discuss connections to recent INS experiments as well as 
other applications of our novel approach. 

Model. We consider a Heisenberg magnet of coupled 
pairs of spins 1/2, with the general Hamiltonian 



(1) 



{ii')r 



defined on a lattice of dimer sites i. J and K are the 
intradimer and interdimer couplings, and m = 1,2 la- 
bels the two spins of each dimer. Without quenched dis- 
order in the couplings J and K, this Hamiltonian can 
describe, e.g., spin ladders and bilayer Heisenber g mag- 
nets, but also applies to materials like TlCuCls 2l|-24 
and BaCuSiaOe U^Mi- 

Generalized bond operators and inhomogeneous con- 
densates. We denote the four states of each dimer i 
by \tk)^, fc = 0,...,3, where |to) = (| U) - I it))/V2, 
Ih) = (-1 n> + I U))/A 1^2) = *(| tt)+ U))/V2, 
Its) = (I tl) + I lt))/^/2■ Formally, bosonic opera- 
tors tlf, can be introduced which create these states out 
of a fictituous vacuum, \tk)i = t|j,|uac)i, with the con- 
straint J2k Ak^ik = 1 defining the physical Hilbert space 



[16j. In a paramagnetic phase, it is convenient |19IJ to 
fully "condense" the singlet, such that the operators 
{a ~ 1,2,3) now create triplet excitations ("triplons") 
from a singlet background, and the Hilbert-space con- 
straint becomes of hard-core type, X^a^a^ia < 1 
Approximations can now be understood as expansion 
about the singlet product state ji/'o) = Di l^o)i- 

Magnetically ordered phases correspond to a conden- 
sate of triplets, such that the reference state now involves 
a linear combination of singlet and triplets on each site. 
A consistent description of excitations requires a basis 
rotation in the four-dimensional Hilbert space spanned 
by the \tk) [H, [2^. Here we generalize the approach of 
Ref. to inhomogeneous states. For every dimer site, 
we introduce a SU(4) rotation to new basis states and 
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operators, 



(fc,fc' = 0,...,3) (2) 



such that iV'o) — Yii l^o)i replaces \tpo) as the reference 
state. For instance, a Nccl state polarized along z is 
obtained from |fo) = (|to) + |t3))/y2 = | U)- 

To make progress, the are chosen such that \ipo) is 
the best product-state approximation to the ground state 
of ■ This can be efficiently found using dimcr mean-field 
theory, by solving the coupled set of mean-field problems 
for each dimcr, 



(S.'™'), (3) 



where the expectation value {Si'm') is taken with |to)i'- 
This is equivalent to a saddle-point approximation, for 
details see Ref. H^. 

To determine magnetic excitations, one now re-writes 
the Hamiltonian ([T|) in terms of the \tk)i using the trans- 
formation ([2]) [s^l. In analogy to the paramagnetic case, 
one fully condenses <J, such that the ij, (a = 1, 2, 3) de- 
scribe excitations on top of the reference state iV'o)- The 
Hamiltonian takes the formH = Hq+Hi+H2+'H3+H4, 
where Tin contains n ia"^ {a = 1,2,3) operators, and 
an additional hard-core constraint for the t,^ is implied. 
With the proper (saddle-point) choice of the reference 
state. Hi vanishes, and ^2 describes Gaussian fluctua- 
tions around a (inhomogeneous) magnetic state. 7^2 has 
the form 



(4) 



where quenched disorder enters via random Ai 

T-L2 is solved by a bosonic Bogoliubov transformation, 
which yields the 3A^ positive-energy eigenmodes of 7^2, 
with N the number of dimer sites. In both paramag- 
netic and coUinearly ordered phases, the three polariza- 
tion directions decouple, giving rise to triply degenerate 
modes in the paramagnetic case and a longitudinal and 
two transverse modes in the collincar case. The eigen- 
states of 7^2 are used to calculate the T = dynamic 
spin susceptbility, as described in Ref. [S^. 

The excitations described by this method interpolate 
continuously between triplons of a paramagnet and spin 



waves of a semiclassically ordered state [281 . I31| . An- 
harmonic effects are neglected at this level of approxi- 
mation; near criticality, this is qualitatively permissible 
if the anomalous exponent 77 is small [77 = 0.03 for the 
Heisenberg model in (2 -|- 1) dimensions]. 

Bilayer magnet. In the remainder of the paper, we 
illustrate the application the method to a simple case, 
namely the bilayer Heisenberg model, as realized e.g. in 
BaCuSi2 06 [25-27|. Here, the dimcrs live on a 2d square 
lattice with interlayer coupling J and intralayer coupling 



K = = Kf^ for nearest- neighbor sites . In the ab- 
sence of disorder, this model is in singlet ground state for 
J ^ K and an AF ground state with ordering wavevec- 
tor Q = (7r,7r) for J <^ K, with a 0(3) quantum critical 
point at iJ/K)c = 2.5220(1) In the harmonic bond- 
operator approach, this transition occurs at {J/K)c = 4; 
including interactions between the triplons allows one to 
obtain a value very close to the QMC result [l^ . 

With disorder of random-mass type, the character of 
the QPT changes, as dictated by the Harris criterion [3^ . 
Numerical simulations have shown that a new critical 
point with conventional power-law singularities obtains 
[33t for not too strong disorder [33| • Here, we shall focus 
on the excitation spectrum in the vicinity of this criti- 
cal point upon inclusion of disorder where signatures of 
quantum Griffiths behavior can be expected 0] . 

Modelling disorder. We shall mainly employ bimodal 
distributions of coupling constants, being experimentally 
relevant to cases where chemical substitution modifies 
exchange paths, as occurs, e.g^ in Tli_a;Ka:CuCl3 [3^ or 
(C4Hi2N2)Cu2(Cli_:j;Bra,)6 [l||. Provided that all cou- 
plings remain AF, this type of bond randomness does 
not introduce magnetic frustration, such that the mag- 
netic order realized for J <C /\ continues to be collinear 
with wavevector Q = (7r,7r). We will denote the corre- 
sponding order parameter, the staggered magnetization 
per spin, by Mg. 

Results: Weak intradimer bonds and evolution across 
QPT. We now present numerical results of our approach, 
for disordered intradimer coupling J, with values Ji and 
J2 taken with probabilities (1— p) and p, respectively, and 
use the interdimer coupling K as a, tuning parameter to 
access the transition. Energies will be quoted in units of 
Ji , the system size is N = L'^ with L = 64 and additional 
2^ supercells (soj . disorder averages are performed over 
~ 50 realizations, and the temperature is T = unless 
noted otherwise. 

Results for small concentrations p of weak intradimer 
bonds, J2 = ./1/2, can be found in Fig. [T] Panels (a,b) 
display Mg as function of K and p, while panels (c)-(h) 
show the dynamical susceptibility x"(9, w) for different 
parameter sets. As expected, this type of substitution 
drives the system towards the ordered phase, by decreas- 
ing the apparent gap and shifting the QPT to smaller K 
upon increasing p. There is a broad range of parameters 
with weak magnetic order, i.e., small non-zero Af,, to be 
discussed below. 

Introducing a small amount p of disorder into the para- 
magnetic phase of the clean system causes a transfer of 
spectral weight to low energies. Fig. [I{c,f-h) and Fig. [5] 
Doping creates excitations inside the gap which eventu- 
ally induce weak order upon increasing p. Importantly, 
this disorder-induced low-energy response, although cen- 
tered around the ordering wavevector Q, is broad in both 
momentum and energy (except for Goldstonc modes at 
extremely small w): This reflects the localization ten- 
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FIG. 1: (color online) Numerical results for the disordered bilayer Heisenberg model, with a concentration of p weak interlayer 
couplings with J2 = Ji/2. a) Phase diagram, showing the staggered magnetization Ms as function of K/Kc for different levels 
of disorder p, with Kc = Ji/4 denoting the critical coupling of the clean system in the linearized bond-operator approach. A 
regime with weak and strongly inhomogeneous order emerges for small p and K < Kc- b) Ms as function of doping level p 
for different K. c-h) Transverse dynamic susceptibihty x" (?*, '-^) along a path in the 2d Brillouin zone [q^ = tt) for different 
combinations of K/ Kc and p. The green dashed lines indicate the dispersion of the clean system with the respective K. 



dency of the corresponding modes, see Fig. [31 For in- 
termediate p near 1/2, the spectrum visibly separates 
into upper and lower branches which exist over the en- 
tire Brillouin zone. Fig. [TJe); these branches correspond 
to modes which are primarily carried by either the Ji or 
the J2 bonds. 

We now turn to a more detailed analysis of the weak 
substitution-induced magnetic order. Fig. [3] portraits a 
single realization of disorder, by showing the spatial dis- 
tribution of J and the local magnetization Mi as well as 
the wavefunctions for selected eigenmodes of 7^2 • Several 
features are apparent: (i) The system develops islands of 
non-zero staggered magnetization. Fig. [3l[b) , with a char- 
acteristic length scale of ^, the correlation length of the 
clean reference system. These islands, being rare events 
in the sense of Griffiths 0], exist in regions with larger 
concentration of weak interlayer bonds, (ii) The distribu- 
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FIG. 2: (color online) Momentum-integrated transverse sus- 
ceptibility x"(^)j for the bilayer model with random weak 
interlayer couplings as in Fig. [1] Disorder tends to close the 
spin gap by transferring of spectral weight to low energies. 



tion of local magnetization values becomes broad on log- 
arithmic scales, Fig.[3l^f), a typical fingerprint of Griffiths 
behavior, (iii) The low-energy "magnon" modes appear 
strongly localized on individual magnetization islands, 
Fig.[3l^c,d). This behavior persists for essentially all ener- 
gies below the gap of the clean reference system, whereas 
higher-energy modes tend to be delocalized. Fig. EJe). 
We have confirmed this localization tendency by analyz- 
ing the inverse participation ratio (30l |. 

Taken together, this implies an interesting "dual" na- 
ture of the weakly ordered phase: It displays a sharp 
Bragg peak in the static structure factor S{q) [s^ at 
Q = (tt, tt), as the underlying order is perfectly staggered 
(albeit strongly inhomogeneous). At the same time, the 
small-w dynamic response w) is anomalously broad 
in g-space, due to the disorder-induced mode localiza- 
tion. Our numerical data support this picture: Fig. IH^a) 
displays S{q) and its finite-size scaling behavior, confirm- 
ing the existence of a sharp Bragg peak, while Fig. WO^) 
shows x"(g, w) at fixed small w, where a spin- wave-like 
two-peak structure is only visible for w/J < 0.02. 

We note that a somewhat similar localization of low- 
energy models has been observed for spin waves in ran- 



domly depleted 2d antiferromagnets |36l . I37IJ ; one impor- 
tant difference to our problem is the presence of random 
Berry phases in the site-depleted antiferromagnet which 
renders disorder stronger than in our random-mass case. 

Results: Other types of randomness. We have studied 
various other types of disorder, with selected results to 
be found in the supplemental material [s^l- Depending 
on the type of disorder, the overall triplon bandwidth 
may either increase or decrease with doping; the same 
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FIG. 3: (color online) Bilayer Heisenberg magnet with 
K/Kc = 0.92 and a single disorder realization of p = 0.02 
weak intradimer bonds with J2 = Ji/2, as in Fig. [UJg), with 
L — 64. (a) Spatial distribution of J values, indicating the 
location of the weak bonds, (b) Corresponding spatial dis- 
tribution of the local magnetization Mi. (c)-(e) Wavefunc- 
tion amplitudes of selected eigenmodes of 'H2. (f) Disorder- 
averaged probability distribution of local (staggered) magne- 
tization values for parameters as in (a)-(e). 



applies to the apparent gap in x"{^)- Band splitting like 
in Fig. [Ije) is visible for sizeable bimodal disorder; for 
continuous distributions this is replaced by strong smear- 
ing of intensity. Most importantly, spectral broadening 
at low energies appears generic; this effect is strongest if 
the dopants tend to drive the system towards the ordered 
phase, as is the case in Figs. [T] and S3 [soj . 

Comparison to experimental results. A few recent ex- 
periments 1^ I13I llSj l have studied bond randomness by 
ligand doping in spin-ladder materials. INS shows exci- 



tations with a larger energy width and an increased spin 
gap as compared to the undoped case. This consistent 
with our results for dopants which create locally stronger 
intradimer couplings. Fig. S2, or weaker interdimer cou- 
plings. Fig. S4, or both. For the investigated materials, 
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FIG. 4: (color online) (a) Static structure factor S{q) and (b) 
susceptibility x\q,u}) at different uj, for K/Kc ~ 0.92 and 
p = 0.02 weak intradimer bonds for L = 128. The inset in (a) 
shows the finite-size scahng of [S{Q)/{2N)]^^'^ which equals 
the order parameter Ms in the thermodynamic limit. 



such properties can indeed be deduced from studies of 
the end members of the doping series (soj . 

While INS data of bond-disordered nearly-critical 
dimer magnets in 2d or 3d are not available to our 
knowledge, it is interesting to connect our results to 
La2-a;Sra;Cu04: At a; = 0.145 this compound displays 
a field-driven magnetic ordering transition, where INS 
data indicate the closing of a spin gap, but in the pres- 
ence of a non-diverging correlation length [ll[. Since it 
is known that La2-2:Sr2;Cu04 shows spatially disordered 
charge stripes [sst which in turn lead to modulated mag- 
netic couplings, our modelling of random-mass disorder 
can be expected to qualitatively describe the soft-mode 
behavior of La2-a;Sr2;Cu04. Therefore we propose that 
quenched disorder tends to localize low-energy magnetic 
modes at the QPT in La2-2;Sra;Cu04, thus providing a 
cutoff to the apparent magnetic correlation length. Dis- 
order is also expected to cut-off one-parameter scaling in 
x"(w,T) - this is testable in future experiments. 

Summary. For coupled-dimcr magnets, we have pro- 
posed an efficient method to calculate real-frequency ex- 
citation spectra in the presence bond randomness across 
the entire phase diagram. This allows one in particu- 
lar to study magnetic excitations near quantum criti- 
cality, where the effect of disorder is generically strong. 
Using the bilayer Heisenberg model as an example, we 
have studied disorder-induced spectral weight transfer 
and weak magnetic order. We find strong broadening 
of low-energy excitation spectra due the localization ten- 
dency of the relevant modes. Further high-resolution INS 
experiments, e.g., on Tli_j;K^CuCl3 |35|, which could 
test our predictions are called for. 

Our method can be extended to the case with mag- 
netic field, in order to access excitations of disordered 
magnon Bose condensates, Bose glasses, or disordered 
supersolids. A generalization to frustrated dimer lattices 
and to incommensurate order is possible as well. 
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I. BOND OPERATORS WITH DISORDER: 
NUMERICS 

In the foUowing wc describe details of our bond- 
operator approach to dimer magnets with quenched dis- 
order, in particular the numerical implementation. We 
work with finite-size systems of = dimers on a 
square lattice with periodic boundary conditions; the 
generalization to other lattices is straightforward. 



A. Optimum product state 

The basis rotation in the four-dimensional Hilbert 
space of each dimer, Eq. (2) in the main text, is deter- 
mined by the requirement that \iPq) = |to)i is the 
best product-state approximation to the ground state 
of %. The state jV'o) can be obtained in various ways: 
(i) by minimizing "Ho = (■(/'o|^|'0o), (ii) by demanding 
■Hi to vanish, as the latter would create a condensate 
of ta bosons, and (iii) dimer mean- field theory, Eq. (3). 
It is easy to see that these methods are equivalent, be- 
cause upon minimizing "Ho one finds a saddle point where 
by construction linear variations (and thus T-Li) vanish. 
Finally, the minimum condition can be converted into 
mean-field theory: The variation of ('0o|^|V'o) w.r.t. \tf))i 



(i) 

or, equivalently w.r.t. Uq^!, reads 



su, 



(i) 



J, 



Ok' 



su, 



Ok' 



E 



su, 



Ofc' 



(SI) 



where (. . .)i is a shorthand for i(io| ■ ■ • \to)i and him = 
J2^'m' K:^''"'{Si'm')^'■ Thc Condition Eq. ^ is the same 
as the minimum condition for the mean-field Hamilto- 
nian 

•H^'f in Eq. (3) of the main text; for the latter 
single-dimer problem the variation of course yields the 
exact ground state. This establishes the equivalence of 
the three criteria. We also note that jV'o) becomes the ex- 
act ground state of Ji in the limit of infinite coordination 
number between the dimers. 

To practically determine the optimum collinearly or- 
dered product state for the Hamiltonian Ji with a given 
realization of disorder, i.e. given random couplings, we 
have found it most efficient to iterate the mean-field prob- 
lem Eq. (3) until convergence. We start with an initial 
collinear guess for the {Si„i). Then, diagonalizing the 
4x4 Hamiltonian H^p in Eq. (3) for each dimer site i 
gives a 4 X 4 matrix of eigenvectors which we choose to 



be [/^*^ In particular, the lowest-energy eigenvector Uq]^, 
can be used to calculate new T = values of the (Sim), 
which closes the self-consistency loop. 



B. Conversion of the Hamiltonian 

To calculate the excitation spectrum, the Hamiltonian 
H needs to be expressed in terms of the local tia excita- 
tion operators, as described in the following. 

The starting point is the representation of the spins 
Sim via transition operators between the states \tk)i of a 
dimer. 



E 

kk' 



(S2) 



with 4x4 matrices s"™ for the spin components 5" (a 
x,y, z = 1, 2, 3) of the m = 1, 2 spins: 
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This is of course equivalent to the bond-operator repre- 
sentation of Sachdev and Bhatti 



ieapytyt^-y) , (S4) 



in the Hilbert space with J^k ^Ik^ik = 1 bosons. After 
the basis rotation (2), Eq. (jS2p becomes 



^im — ''^^f,kk'\^k)i i{ik' 
kk' 



(S5) 



with the transformed spin matrices now being site- 
dependent: 



^i,kk' 



(S6) 



2 



The full Hamiltonian in terms of transition operators 
takes the (exact) form 

T~L = X^^ \ik)i i{ik'\ 

ikk' 

+ E Y^y'''\ik)^\ii), (S7) 



ii'kk'll 
^kk' 



Its coefficients X^'' and Fjf^ " are related to the original 



J 



couplings Ji and -K'™'" , respectively. They can be gener- 
ated numerically for arbitrary u[2> of a finite-size system 
with a given realization of disorder, using the equations 
(1), jSSl), and ([S6|). 



This Hamiltonian is now rewritten into Ua {a = 1,2,3) 
bosonic operators which are defined to act on the back- 
ground state nj^o)i, i-e., iLl^o)* = \ia)i- Using the 
decomposition H = Ho + Hi + T-L2 + T-Ls + "Ha we have 



Ho 
Hi 

H2 



i ii' 
ial3 ii'ap 

I 



ii' ol(5 



(S8) 



where a, (3 ~ 1,2,3. Similar expressions can be found 
for H3 and H4 which are, however, not needed in the 
following. In each run, we have checked that Hi vanishes 
to numerical accuracy provided that the iteration of the 
product state has converged. 



C. Gaussian fluctuations and susceptibility 

The Gaussian fluctuations of interest are contained in 
H2, which describes 3N modes where N is the number 
of dimers, and the factor of three arises from the spin 
polarizations a = x,y, z. 

We generate a matrix form of H2 numerically from the 
coefficients and Y-f^^' via Eq. In both the 

paramagnetic and the coUinearly ordered cases studied 
here, the polarization directions decouple, such that H2 
is block-diagonal. Then, it remains to perform a bosonic 
Bogoliubov transformation for the transverse and lon- 
gitudinal fluctuations separately. Each amounts to the 
numerical diagonalization of a 2N x 2N non-Hermitean 
matrix. We have implemented the procedure described 
in the appendix of Ref. H using standard LAPACK rou- 
tines. An efficient implementation can be achieved using 
the OpenMP version of Intel's MKL library, which al- 



128'^ 



the 



lowed us to access system sizes up to N 
coUinear case using a desktop computer. 

Each diagonalization yields (per polarization) 2A^ 
eigenvectors and eigenvalues, with the latter coming in 
pairs with opposite sign. In the presence of a symmetry- 
breaking condensate, we always find at least one pair of 
eigenvalues with magnitude smaller than 10~^J in the 



transverse piece of H2, corresponding to the expected 
zero-energy (Goldstone) mode. 

The positive-energy eigenvalues uj„ and their eigenvec- 
tors (u„, Vn) describe the Gaussian fluctuation modes t„ 
of H, which are linear combinations of the t operators: 



"nia '^iol 



(S9) 



The LUn and (u„,?;„) can be used to extract physical ob- 
servables. In the paper, we have concentrated on the dy- 
namic spin susceptbility at T = in the single-magnon 
approximation 



J2\M;:{q)\'5{u 



(SIO) 



where []av denotes disorder averaging, and M^{q) is the 
matrix element {0\S"{q)\n) of the Fourier-transformed 
spin operator (in the single-magnon approximation). 



with the eigenmode 



r^|0). Here, nm is the lattice 



coordinate of an individual spin. For the bilayer model 
under consideration, (f has a z component Qz = 0,7r; we 
restrict our attention to the staggered response, qz = 
TT, because the low-energy modes have significant weight 
only in this channel. 

For a L X L system, Xa('?*, w) can then be calculated 
with a momentum resolution in {qx,qy) of L points per 



direction. The delta function in Eq. (jSlOp may be con- 
verted into a Lorentzian with artificial width. However, 
for momentum-resolved spectra we found it advantageous 
to use binning instead, because the exact intensity in the 
clean case diverges as l/w upon approaching Q for trans- 
verse (Goldstone) modes, such that a Lorentzian broad- 
ening leads to unphysical weight at elevated energies. 

In order to improve the momentum resolution, we have 
also implemented twisted boundary conditions, or equiv- 
alently, supercells: The system is treated as a periodic 
arrangement of Ad x M subsystems of size L x L each 
and overall periodic boundary conditions. Then, di- 
agonalizations of a system with size are required, and 
the number of momentum points per direction is M x L. 
We note that supcrccU calculations render the Hamilto- 
nian matrix inevitably complex, while it is otherwise real 
in the coUinear case considered here. Most results in the 
paper have been obtained for L = 64 and M = 2, and we 
show the transverse piece of x" only. 



D. Static structure factor 

The static spin structure factor is defined as 



Siq) 



1 

m 



(S12) 



with 2N the total number of spins and fu'mm' = '^im ~ 
film' the distance between a pair of spins. In the bi- 
layer model, the ordered state displays a Bragg peak at 
Q = (7r,7r,7r); for iV — > oo its intensity is related to the 
staggered magnetization per spin via = S{Q)/{2N). 

In Fig. 4 of the main text we show S{q) for Qz = tt; for 
the finite-size scaling S{q) has been calculated in a mean- 
field approximation only, {Sim ■ Si>m') -> {Sim) ■ {Si'm')- 
For small systems we have checked that the fluctuation 
corrections only lead to minor quantitative changes. 



II. BOND OPERATORS FOR THE CLEAN 
BILAYER MAGNET 

To make the presentation self-contained, we quickly 
summarize the most relevant results of the bond-operator 
approach to the bilayer Heisenberg magnet, Fig. lSli with- 
out quenched disorder. 



A. Non-interacting bosons 

A number of papers have applied versions of the bond- 
operator formalism to the bilayer Heisenberg model. In 
many cases, interactions between triplet excitations are 
neglected, i.e., only bilinear terms in the bond-operator 
Hamiltonian are kept. The approaches differ in the treat- 
ment of the Hilbert-space constraint: It is either fulfilled 




FIG. SI: Bilayer Heisenberg model with interlayer coupling 
J and intralayer coupling K. 



on average, leading to so-called bond-operator mean-field 
theory,—^— or ignored completely (at the quadratic level) 
in the spirit of linear spin-wave theoryj^i^ Here we focus 
on Rcf. Is which follows the latter concept, as this allows 
one to consistently treat the ordered phase across the 
entire phase diagram such that Goldstone's theorem is 
automatically fulfilled. We also restrict ourselves to the 
case of zero applied field. 

The basis transformation in Eq. (2) in the main text 
takes the following simple form: 



''iO 



4 + Ae''3-«.4 



'-a 



andili =t\i.tl2 



1 + A2 
rl 



'•13 



TTa^ (S13) 



t\-^ . Here A is parameterizes the linear 
combination of singlet and z triplet, and Q = (7r,7r) the 
ordering wavevector in the dimer lattice of the Ri . Min- 
imization of the product-state energy, Hq = (■0o|^|V'o)i 
yields 



K -U 
K + AJ' 



(S14) 



For K/J < 1/4 the system is paramagnetic, and from 
'H2 one finds the energy dispersion to be 



7g- = (cosq.^ + cosqj^)/2, 
with a spin gap given by 



(S15) 
(S16) 

(S17) 
(S18) 



For K/J > 1/4 the system is in an antifcrromagnetic 
state, with a product-state staggered magnetization 



A = ^/JiJ -AK). 
and a triplon bandwidth of 



2A 



1-HA2 



(S19) 



per dimer. This value is somewhat reduced by Gaussian 
fiuctuations; in particular in the limit of decoupled layers, 
J — > 0, where A = 1 one recovers the spin-wave value 
Ms = 0.66 (Ref.li). 
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B. Beyond non-interacting bosons 

Going beyond the non-interacting (i.e. linearized) the- 
ory requires to take into account the quartic interactions 
in ^4 (cubic terms arc absent from the bilayer model) as 
well as the hard-core constraint. An efficient approxima- 
tion for the latter is given by the Brueckner theory and 
yields an accurate value for the location of the critical 
point, as shown by Kotov et al.^ However, this approach 
cannot be easily generalized to the ordered phase. An 
earlier approach by Chubukov and Morr— , inspired by 
non-linear spin-wave theory, works in both phases, but 
suffers from divergencies at higher orders, probably be- 
cause it lacks a small control parameter. Thus, we are 
not aware of a tractable and consistent way to incorpo- 
rate interactions into the bond-operator approach for all 
phases. 



III. RESULTS FOR DIFFERENT TYPES OF 
DISORDER 

In the main paper, results were discussed for disor- 
dered interlayer (or intra-dimer) couplings J, specifically 
a small concentration of weakened interlayer couplings. 
Here we first discuss some general aspects of disorder, 
also with an eye towards real materials, and then present 
data for different types of disorder. 



A. Griffiths regime and paramagnetic gap 

For the bounded disorder distributions considered in 
this paper, there is always a well-defined (hard) exci- 
tation gap for sufficiently small K/ J. Upon increasing 
the gap closes at the boundary to the quantum Grif- 
fiths regimei^ In our mean-field-based treatment, (weak) 
magnetic order sets in at this point; non-linear fluctua- 
tion corrections would be required to obtained the correct 
physical behavior, namely a magnetic phase transition 
inside the Griffiths regime. In other words, at the level 
of our approximation a paramagnetic quantum Griffiths 
regime is absent, and the entire Griffiths regime is mag- 
netically ordered. Despite this, we can expect that the 
properties of the excitations are captured correctly ex- 
cept at ultra-low energies. 

Depending on the type of disorder, the "apparent" en- 
ergy gap in x" may be much larger than the actual gap, 
namely if the intensity in x" is very small for small uj - 
this happens if the low-energy excitations are carried by 
rare configurations (i.e. rare clusters of small J or large 
K). This effect can be seen in Fig. 2a of the main paper, 
where a true gap is absent, but the p = 0.05 susceptibility 
has an apparent gap of sa O.IJ. 



B. Disorder in experiments 

Apart from fundamental interest, specific mod- 
els of disorder may be motivated by the ex- 
perimental situation. Here, bond disorder is 
introduced via ligand substitution, as e.g. in 
Tli.^K.CuCla^'ii (C4Hi2N2)Cu2(Cli_,Br,)6fi^^ 
IPA-Cu(Cl,Bri_,)3^^ or (Hpip)2CuCl4(i_,)Br4,^ 
such that bimodal disorder models appear appropriate. 

In general, such substitutions influence both intra- 
dimer (J) and inter-dimer (K) couplings; this can be 
deduced from the end members of the substitution se- 
ries. For instance, TlCuCls and KCuCla have drasti- 
cally different spins gaps of 0.7 and 2.5 meV, respec- 
tively, and flts to INS data results in values of J = 
5.42 mcV and ^^1,2,3 = (-0.47, -1.43, 0.62) meV for 
TlCuCls, in contrast to J = 4.29 mcV and Ki,2,3 = 
(-0.21, -0.34, 0.37)meV for KCuCls; here i^i,2,3 are 
the strongest couplings in the three-dimensional dimer 
network of these materialsf^i Therefore, a substitu- 
tion Tl— >-K in TlCuCla tends to increase the spin gap; 
the same apphes to Cl^^Br in (C4Hi2N2)Cu2Cl6 and 
(Hpip)2CuCl4. 

It can be expected that partial substitution of ligands 
locally changes J and/or K couplings, but this happens 
in a correlated fashion: A single substituted ligand atom 
influences numerous local bond lengths and angles, such 
that various couplings in the vicinity of the substituent 
are modifled - this applies in particular to interdimer 
couplings which typically result from multiple exchange 
paths. Therefore, models of uncorrelated bimodal disor- 
der are a strong simpliflcation and cannot fully capture 
the disorder physics of real materials. We have therefore, 
in addition to bimodal uncorrelated disorder in J, also 
considered correlated disorder in K, see below. 



C. Bimodal bond disorder 

The case of disordered interlayer couplings J was ex- 
tensively discussed in the main text. There the focus 
was on a small concentration of weakened interlayer cou- 
plings. Those have the tendency to reduce the spin gap 
by inducing in-gap states, which are strongly localized 
and lead to broad modes at lowest energies near the tran- 
sition (Fig. 1). 

In Fig. [52] we show results for the opposite case of a 
small concentration p of stronger interlayer (J) couplings. 
Consequently, the finite-p excitation spectrum now shows 
a larger apparent gap as compared to the p ~ case, 
and ordering tendencies are suppressed. The few strong 
bonds induce non-dispersing excitations above the main 
band and with weak intensity; at the same time the over- 
all width of the main band is reduced, because coherent 
triplon hopping is preempted through the strong bonds. 
Near the transition. Fig. IS2f c') there is still significant 
mode broadening at low energies. 
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(a) KIK, = 0.96, p=0. 1 (b) KIK^ = 1 .0, p=0. 1 (c) KIK^ = 1 .02, p=0. 1 (d) K/K, = 1 .04, p=0. 1 




(0,0) (Tt.Tl) (71,0) (0,0) (7t,7l) (71,0) (0,0) (7l,7t) (71,0) (0,0) (7I,7l) (71,0) 



FIG. S2: Numerical results for the susceptibility x"{<1j^) in ^ bilayer Heisenberg model with disordered interlayer couplings as 
in Fig. 1 of the main paper, but here for density of p strong interlayer couplings with J2 — 3Ji/2. 




(0.0) (7C,7l) (71,0) (0,0) (71,71) (7t,0)(0,0) (7l,7t) (71,0) (0,0) (7t,7l) (71,0) 



FIG. S3: Susceptibility x"(<?i but now for disorder in the intralayer couplings K: a density of p dimers has strong intralayer 
couplings K2 = 2A'i with all neighbors. 




(0,0) (7t,7l) (71,0) (0,0) (71,71) (71,0) (0,0) (7I,7t) (71,0) (0,0) (7t,7l) (7t,0) 



FIG. S4: Susceptibility x"{q,Lo) for disordered intralayer couplings K as in Fig. IS3I but now a density of p dimers has weak 
intralayer couplings with K2 ~ A'i/2 with all neighbors. 



We now turn to disorder in the intralayer (or intcr- 
dimer) couplings K. For uncorrelated disorder in K, we 
have found that disorder effects are much weaker than 
that of uncorrelated disorder in J, i.e., uncorrelated K 
disorder efficiently averages out because of the presence of 
many (here eight) K bonds for each dimer. Considering 
that uncorrelated disorder is unrealistic for real materi- 
als, we have chosen to model correlated disorder in the 
following fashion: We randomly choose a (small) frac- 
tion p of dimers, mimicking the location of dopants, and 
then give all K bonds emerging from such dimers the 
strength K^j while the remaining K bonds have strength 
Ki. Results are shown in Figs. IS3I and IS4I for the cases of 
p dimers with stronger and weaker intralayer couplings, 
respectively. 

Clearly, a small concentration of sites with stronger 
intralayer couplings, Fig. IS31 has an effect similar to a 
small concentration of weaker interlayer couplings, Fig. 1: 
The dopants introducde significant in-gap spectral weight 



and drives the system towards the ordered state. The 
in-gap states themselves now show significant energetic 
structure, arising from cluster of two, three, and more 
dimers sites with stronger intradimer coupling. 

In contrast, a small concentration of sites with weaker 
intralayer couplings, Fig. IS41 leads to spectra similar to 
the case of a small concentration of stronger interlayer 
couplings, Fig. IS21 with bandwidth reduction and mode 
broadening being the dominant effects. 



D. Box bond disorder 

For completeness, we have also considered uncorrelated 
continuous distributions of coupling constants. Specifi- 
cally, wc have studied box distributions for either J (with 
non-disordered K) or K (with non-disordered J). Cor- 
responding results are shown in Figs. [S5]and lS61 respcc- 
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(a) ^//f, = 0.88 (h)K/K^ = 0.92 ic)K/K^ = 0.94 (d)K/K^ = 0.96 




(0,0) (Jl,7t) (71,0) (0,0) (71,71) (71,0) (0,0) (71,71) (71,0) (0,0) (7t,7l) (JC,0) 



FIG. S5: Susceptibility x"{Qj^)j but now for a box distribution of interlayer couplings J with a width SJ = J/2, i.e., the Ji in 
Eq. (1) vary in the interval [0.75 J, 1.25 J]. 



ia)K/K^ = OM {b)K/K^ = 0.92 (c)KIK^ = Q.9A {A)KIK^ = Q.96 




(0,0) (71,71) (71,0) (0,0) (7t,7l) (71,0) (0,0) (71,71) (71,0) (0,0) (7t,7l) (71,0) 



FIG. S6: Susceptibility x"('?,'^), but now for a box distribution of intralayer couplings K with a width &K = K, i.e., the K]^, 
in Eq. (1) vary in the interval [Q.bK, 1.5A']. 



tively. 

Not surprisingly, box distributions tend to enlarge the 
effective bandwidth of magnetic excitations, simply be- 
cause the range of available couplings increases: For in- 
stance, a variation of J tends to shift the center energy of 
the triplon band, such that a system with box-distributed 
J . Fig. IS5[ effectively samples bands corresponding to a 
range of J values. Similarly, a variation of K tends to 
modify the bandwidth, such that the bandwidth for box- 
distributed K G [K-5K/2, K+5K/2] is larger than that 
corresponding to the pure-ii' case. However, the increase 
in bandwidth is moderate even for large 5K, Fig. IS61 
which demonstrates that uncorrelated disorder in the in- 
tradimcr couplings averages out as noted above. 

Smearing of the excitations (but no splitting) is visible 
at all energies for both box distributions, with smearing 
again being most pronounced at low energies near the 
transition. In particular, strongly broadened low-energy 
modes are visible in Figs. [S5fc and[ 



IV. 



EVOLUTION OF MODE BROADENING 
AND BAND WIDTH 



The calculated x"(9i ^) enables a quantitative analysis 
of the evolution of spectral properties with the disorder 
level, in particular the energetic mode broadening and 
the triplon bandwidth. 

A first impression can be obtained by considering 
the energy-dependent response at fixed wavevector, as 
shown in Fig. ISTf a-b) for the putative ordering wavevec- 



tor q = Q= (tt, tt). Both panels are for different disorder 
in the interlayer couplings: (a) has a small concentration 
of strong interlayer bonds, which have the simple ten- 
dency to broaden the mode and shift it to higher energies. 
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FIG. S7: (a) Excitation spectrum ■x'{q,ui) at fixed q = Q = 
(7r,7r), for a small concentration of strong interlayer bonds 
with J2 = 3J1/2, as in Fig. IS2I (b) Same as (a), but for a 
small concentration of weak interlayer bonds with J2 = Ji/2, 
as in Fig. 1 of the main paper. (c,d) Linewidth and triplon 
bandwidth for the case of of strong interlayer bonds with J2 = 
3J1/2 of function of p. The error bars arise primarily from 
the limited energy resolution of our finite-size calculation; the 
lines are guide to the eye. 
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In contrast, (b) has a small concentration of weak inter- 
layer bonds, which have two effects: The "main" mode 
is slightly shifted upwards, but a significant amount of 
weight is created inside the gap. 

In the case (a) of few strong interlayer bonds, it makes 
sense to quantify the mode broadening, this is done in 
Fig. IS7f c) for two different values of K/Kc correspond- 
ing to the data depicted in Fig. [S2ja,b). The data are 
consistent with a linewidth scaling linearly with p, with 
a prefactor which increases upon approaching the QPT. 
(As noted before, the broadening is largest in the regime 
where the apparent spin gap closes, see Fig. ISSlc).) 

Fig- ISTT d) displays the evolution of the magnon band- 
width for the same case of few strong interlayer bonds; 
the small-intensity non-dispersing states at the upper end 
of the spectrum, Fig. lS21 have been ignored in this analy- 
sis. The bandwidth decreases as expected from Eq. (jSlSP : 
the decrease is approximately linear for small p. This is 
qualitatively matches recent experimental INS dataJ^ii^ 

We finally note that non-linearities, i.e., triplon-triplon 
interactions, will also lead to broadening and renormal- 
ization effects, which are not part of the present calcu- 
lation. However, interaction-induced broadening is typ- 
ically small at low energies, i.e., near the bottom of the 
excitation band, because of phase-space restrictions for 
particle decay. Therefore, the disorder effects captured 
by our calculations can be expected to dominate at low 
energies. 



V. LOCALIZATION OF LOW-ENERGY MODES 

As illustrated in Fig. 3(c,d) of the main text, the 
low-energy excitation modes of the disordered coupled- 
dimer magnet tend to be spatially localized. This can be 
quantified by considering the inverse participation ratio 
(IPR) of the eigenmodes, with the general definition of 
IPR — for a wavefunction !■(/;), and V'i = {A'^) 

where i is a lattice site. The finite-size behavior of 
the IPR is indicative of localization: For spatially ex- 
tended states, the IPR scales with the system size as 
whereas it approaches a constant for localized states. For 
exponential localization, the IPR value is then a measure 
of the localization length, IPR oc 

For the eigenvectors t„ of the bosonic Bogoliubov 
transformation (IS9I) we define the IPR as 



L = .32 



L = 64 



L = 128 



IPR„ = ^ (|WmQp - \v„ia\'^y 



(S20) 



note that ^ 



l^niap) = 1. Some previous 



papers* have employed an alternative definition. 



ipr; 



(S21) 



we have checked that both quantities show qualitatively 
similar properties. 




1 2 3 4 5 



FIG. S8: Inverse participation ratio (jS20[) for the eigen- 
modes of 14,2 of the disordered bilayer Heisenberg model, 
with KjKc = 0.92 and p = 0.02 weak interlayer bonds with 
J2 = Ji/2 as in Fig. 1(g) of the main paper. The panels show 
the IPR as function of the mode energy for different system 
sizes L X L; the color encodes the average number of modes 
(per IPR and energy interval and per disorder realization) in 
a histogram fashion. Localization tendencies, with the IPR 
being independent of system size, are pronounced for energies 
uj < 0.3J, i.e., below the gap of the clean system. 



L=32 



L = 64 



L = 128 



^ 0.01- 



fyj 




FIG. S9: Same as Fig. IS8I but for a box distribution of in- 
tralayer couplings K with K/Kc = 0.92 as in Fig. [S6|b). 



Energy-resolved histograms of IPR values (|S20|) are 
shown in Fig. lSSI for the same parameters as employed in 
Figs. 1(g), 3, 4 of the main text, i.e., p = 0.02 weak inter- 
layer bonds. A scaling with inverse system size is visible 
for essentially all energies to > 0.35 J, indicating extended 
states, whereas states with lu < 0.3J are clearly local- 
ized. Considering that the spin gap of the clean system 
is A = 0.28 J, it is obvious that it is the disorder-induced 
states below the gap which are subject to localization. 
(The somewhat enhanced IPR at both w/ J « 1 and the 
upper band edge arises from the flat pieces of the mode 
dispersion.) Comparing Figs. 1(g) and [S8l also shows that 
energies with a large IPR correlate with those showing a 
large broadening in momentum space, both indicative of 
a short real-space localization length. 

Fig. [S9] illustrates similar localization tendencies for 
a box distribution of couplings, albeit with a somewhat 
larger localization length. Here, localization takes place 
a both band edges, but concerns a broader energy range 
near the band bottom as compared to the top. 

We conclude that disorder-induced localization of low- 
energy modes is generic near the QPT of the bilayer 
Heisenberg model. 
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